Received June 18, 2010; Accepted July 7, 2010 

Preprint typeset using I^TeX style emulateapj v. 11/10/09 



DIMENSION AS A KEY TO THE NEUTRINO MECHANISM OF CORE-COLLAPSE SUPERNOVA 

EXPLOSIONS 

J. NoRDHAus^ A. Burrows^ A. Almgren^, J. Bell^ 

Received June 18, 2010; Accepted July 7, 2010 

ABSTRACT 

We explore the dependence on spatial dimension of the viability of the neutrino heating mechanism of 
core-collapse supernova explosions. We find that the tendency to explode is a monotonically increasing 
function of dimension, with 3D requiring ~40— 50% lower driving neutrino luminosity than ID and 
~15— 25% lower driving neutrino luminosity than 2D. Moreover, we find that the delay to explosion 
for a given neutrino luminosity is always shorter in 3D than 2D, sometimes by many hundreds of 
milliseconds. The magnitude of this dimensional effect is much larger than the purported magnitude 
of a variety of other effects, such as nuclear burning, inelastic scattering, or general relativity, which 
are sometimes invoked to bridge the gap between the current ambiguous and uncertain theoretical 
situation and the fact of robust supernova explosions. Since real supernovae occur in three dimensions, 
our finding may be an important step towards unraveling one of the most problematic puzzles in 
stellar astrophysics. In addition, even though in 3D we do see pre-explosion instabilities and blast 
asymmetries, unlike the situation in 2D, we do not see an obvious axially-symmetric dipolar shock 
oscillation. Rather, the free energy available to power instabilites seems to be shared by more and 
more degrees of freedom as the dimension increases. Hence, the strong dipolar axisymmetry seen in 
2D and previously identified as a fundamental characteristic of the shock hydrodynamics may not 
survive in 3D as a prominent feature. 

Subject headings: hydrodynamics - supernovae: general - stars: interiors - neutrinos 



1. INTRODUCTION 

It was shown some time ago that multi-dimensional 
instabilities obtain and are probably central to the core- 
collapse supernova mechanism (Burrows & Fryxell 1992; 
Herant, Benz, & Colgate 1992; Herant et al. 1994; Bur- 
rows, Hayes, & Fryxeh 1995; Janka & Miiller 1996). 
However, it has not been clear whether modest effects, 
particularly in the neutrino sector, that individually 
might have little consequence, could accumulate to col- 
lectively push the core beyond the threshold of instabil- 
ity into explosion (Mezzacappa et al. 1998,2001,2007; 
Bruenn et al. 2007,2010; Marek & Janka 2009). More- 
over, though the explosion energy is only 10^^ ergs, since 
more than 10^^ ergs of neutrinos issue from collapse, 
the supernova phenomenon might seem to be a "one- 
percent" effect, necessitating attention to every detail to 
see the qualitative effect that is the supernova. However, 
in the crucial early post-bounce phase of a few hundred 
milliseconds, only a few times 10^^ ergs in electron-type 
neutrinos emerges, making the neutrino-driven explosion 
more like a "tens of percent" effect, challenging the no- 
tion that fine detail in the neutrino sector is the key to 
the viability of the neutrino explosion mechanism. 

It may well prove to be the case that the fundamental 
impediment to progress in supernova theory over the last 
few decades has not been lack of physical detail, but lack 
of access to codes and computers with which to properly 
simulate the collapse phenomenon in 3D. This could ex- 
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plain the agonizingly slow march since the 1960's towards 
demonstrating a robust mechanism of explosion. State- 
of-the-art two-dimensional simulations are still ambigu- 
ous and problematic (^21), though they manifest instabil- 
ity and turbulence anoincrease the dwell time of matter 
in the gain region (Bethe & Wilson 1985), and, hence, 
the efficiency of the neutrino-matter coupling (Murphy 
& Burrows 2008). The gain region is the low-optical (- 
neutrino) depth region behind the shock where neutrino 
heating exceeds neutrino cooling. Moreover, relative to 
the spherical case, 2D effects have been shown to enlarge 
the gain region itself and place more matter in shallower 
climes of the gravitational potential well (from which it 
is easier to launch) . Two-dimensional models also enable 
accretion in one quadrant to power explosion in another 
(Burrows et al. 2007b), something impossible in ID and 
one reason neutrino-driven explosions in ID are gener- 
ally thwarted. In ID, an outward explosion diminishes 
the mass accretion rate onto the residual core that pow- 
ers a good fraction of the driving neutrino luminosity 
during the critical incipient explosive phase. 

However, the difference in the character of 3D turbu- 
lence, with its extra degree of freedom and inverse cas- 
cade to smaller turbule scales than are found in 2D, may 
further increase the time matter spends in the gain re- 
gion, so that the critical condition for explosion (Burrows 
& Goshy 1993; Murphy & Burrows 2008) is more easily 
achieved. Our thesis in embarking upon this project was 
that if going from 2D to 3D improves prospects for the 
neutrino-driven explosion by the same degree as already 
demonstrated by going from ID to 2D (see q2|, then the 
neutrino-driven mechanism would perhaps be shown to 
be not only viable, but robust. In that case, the details 
of neutrino interactions, general relativity, and nuclear 
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physics would be of secondary importance for demon- 
strating the mechanism of explosion, though would still 
be centrally important for determining the actual explo- 
sion energies, the magnitude and systematics of pulsar 
kicks, the residual neutron star masses, and the nucle- 
osynthesis, to name just a few important aspects of the 
core-collapse supernova phenomenon. 

Hence, in this paper we have conducted a parame- 
ter study, similar to that performed by Murphy & Bur- 
rows (2008) (see Q, but directly comparing ID (spheri- 
cal), 2D (axisymmetric), and 3D (Cartesian) simulations, 
while varying the driving neutrino luminosity. To isolate 
the effect of dimension, we ensured that all other aspects 
of the runs (equation of state, progenitor model, and neu- 
trino heating algorithm, etc.) were the same. Note that 
for a given massive-star progenitor, the accretion rate in 
front of the stalled shock is independent of dimension, 
since the initial model is spherically symmetric and the 
post-shock matter is out of sonic contact with the super- 
sonic infalling matter in front of the shock. Therefore, 
the character of infall, the evolution of the accretion rate 
and the mass interior to the shock, and by extension the 
evolution of the neutrino luminosities, are functions more 
of progenitor than of dimension (until explosion). 

In what follows, we repeat the ID and 2D study of 
Murphy & Burrows (2008), but, extending it to 3D, ver- 
ify that the viability of explosion by delayed neutrino 
heating is indeed a monotonically increasing function of 
dimension. In going from 2D to 3D, the enhancements 
in the efficiency of the neutrino-matter coupling in the 
gain region and the decreases in the critical luminosity 
are comparable to those seen in going from ID to 2D. 
In a very real sense, this implies that it is ~50% easier 
to explode in 3D than in ID, a huge difference. Since 
real supernovae occur in three domensions, this conclu- 
sion may be a major step in unraveling one of the most 
recalcitrant puzzles in astrophysics. 

In ^ we review the current status of numerical super- 
nova theory. Then, in ^ we summarize the CASTRO 
code and the simple algorithm by which we incorporate 
neutrino heating and cooling for the purposes of this pa- 
rameter study. We follow in ^4] with a discussion of our 
general results, highlighting the role of dimension alone 
in facilitating explosion. We conclude in ^with a sum- 
mary of what we have learned. 

2. SUMMARY OF CURRENT SIMULATION EFFORTS AND 
RESULTS 

Colgate, Grasberger, & White (1961) pioneered the di- 
rect (hydro dynamic) mechanism of core-collapse super- 
novae, in which the bounce shock, launched when col- 
lapse is reversed, propagated outward without stalling. 
Following this, Colgate & White (1966) proposed the 
neutrino mechanism, wherein the agency of explosion 
was a burst of neutrino heating behind the shock after 
bounce. Neutrinos must be included at the high temper- 
atures and densities achieved in collapse. Arnett (1966) 
and Wilson (1971) refined this model, with the latter 
challenging its efficacy. These simulations were done in 
spherical symmetry and employed less sophisticated nu- 
merical techniques and cruder microphysics than can be 
marshalled today. Importantly, current thinking is that 
the direct hydrodynamic mechanism never works. The 
neutrino burst that occurs at shock breakout of the neu- 



trinosphere and, to a lesser degree, photodissociation of 
nuclei into nucleons at the nascent shock debilitate it 
into accretion. Moreover, and disappointingly, modern 
simulations with the latest nuclear and neutrino physics, 
up-to-date progenitor models, and sophisticated numer- 
ics demonstrate that if stellar core collapse is constrained 
to be spherically-symmetric most cores can not super- 
nova (Burrows, Hayes, and Fryxell 1995; Mezzacappa et 
al. 2001; Rampp and Janka 2000,2002; Liebendorfer et 
al. 2001,2005; Buras et al. 2003; Thompson, Burrows, 
and Pinto 2003). However, Kitaura et al. (2006) do ob- 
tain a spherical neutrino-wind-driven explosion for the 
8.8-M0 progenitor of Nomoto & Hashimoto (1988), after 
a slight post-bounce delay (see also Burrows, Dessart, 
& Livne 2007). Such neutrino-driven winds are generic 
after explosion (Burrows, Hayes, & Fryxell 1995), but 
are too weak to generate a 10^^ -erg explosion. This pro- 
genitor can explode in ID because its envelope is ex- 
tremely rarified, and the wind emerges almost immedi- 
ately, but the associated energy is uncharacteristically 
meager (~2.5 x 10^^ ergs, more than an order of mag- 
nitude smaller than the canonical value). It explodes 
spherically because there is no inhibiting accretion tamp 
and can not be generic. 

In the 1980's, Wilson (1985) suggested the "delayed" 
mechanism of explosion, wherein the stalled shock was 
revived by neutrino heating, but after hundreds of mil- 
liseconds. This mechanism was championed by Bethe & 
Wilson (1985), who introduced the concept of the "gain 
region" behind the shock, where net neutrino energy de- 
position was positive. Though his calculations were done 
in ID, Wilson required "neutron- finger" convection be- 
low the neutrinospheres to boost the driving neutrino 
luminosities. Such instabilities, and the corresponding 
luminosity boosts, have been called into question by 
Bruenn & Dineva (1996) and Dessart et al. (2006a). 

In the 1990's, the available computer power enabled 
detailed numerical explorations beyond spherical sym- 
metry into two dimensions (Burrows and Fryxell 1992; 
Herant, Benz, and Colgate 1992; Herant et al. 1994; Bur- 
rows, Hayes, and Fryxell 1995; Janka and Miiller 1996), 
with some of these 2D calculations leading to delayed 
neutrino-driven explosions. These were the ffist simu- 
lations of collapse and explosion that demonstrated the 
existence and potential centrality of hydrodynamic insta- 
bilities and neutrino-driven convection in the supernova 
phenomenon. Fryer and Warren (2002,2004) and Fryer & 
Young (2007) later performed 31) simulations, but used 
an SPH hydro code that may not have adequately han- 
dled the resultant instabilities and turbulence. Never- 
theless, those were some of the ffist 3D simulations in 
supernova theory that attempted to include the relevant 
physics. 

Mezzacappa et al. (1998) challenged the notion 
that multi-D effects were central to the supernova phe- 
nomenon, suggesting that precision spherical simula- 
tions for which detailed neutrino transport was incor- 
porated would be necessary to demonstrate robust ex- 
plosions. However, this is not the current view, with 
multi-D effects ffist explored in the 1990's now assum- 
ing center stage. More recently, using purely hydrody- 
namic calculations, Blondin, Mezzacappa, & DeMarino 
(2003) showed that the shock itself, even without neu- 
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trino heating and the consequent convection, could be 
unstable, introducing the "standing-accretion-shock in- 
stability" (SASI). This hydrodynamic phenomenon has 
been studied analytically and in detail by Foglizzo et al. 
(2006,2007). One can not easily disentangle neutrino- 
driven convection and the SASI, and it may be that 
neutrino-driven convection is predominant (Fernandez & 
Thompson 2009). The major contribution of the Blondin 
et al. paper was to highlight the need to perform 2D cal- 
culations over the full 180° so as not to suppress the 
interesting dipolar component. Most previous calcula- 
tions had been performed on a 90° wedge. Otherwise, 
the "SASI" had been manifest and observed naturally in 
previous multi-D simulations. 

Buras et al. (2006ab), using their code MuDBaTH, ob- 
tained an explosion in 2D with the 11.2-M0 progenitor 
of Woosley & Heger (2006). However, this model's explo- 
sion was underpowered by an order of magnitude. More 
recently, Marek & Janka (2009) obtained an explosion of 
a rotating IS-M© progenitor (their model M15LS-rot), 
a more representative star. However, they experienced 
a long delay to explosion of more than 600 milliseconds, 
by which time the mass in the gain region was too low 
to absorb enough neutrino energy to achieve more than 
~2.5 X 10^^ ergs by the end of their simulation. More- 
over, the shock wave was not followed beyond ^600 km 
and they used a soft nuclear equation of state with an 
incompressibility at nuclear densities, of 180 MeV. 
Their model with K = 263 MeV did not explode. An in- 
compressibility nearer 240±20 MeV is currently preferred 
by measurement (Shlomo, Kolomietz, & Colo 2006). In 
addition, Marek and Janka's high-resolution run of the 
same model (M15LS-rot-hr) did not explode, though it 
was continued to nearly the same final post-bounce epoch 
as their exploding model M15LS-rot. 

The Oak Ridge supernova group, led by A. Mezza- 
cappa and S. Bruenn, is paralleling in many ways the 
work of the Marek & Janka (2009), but is using its 
CHIMERA code (Bruenn et al. 2007,2010; Mezzacappa 
et al. 2007) and Bruenn's ID multi-group, flux-limited 
diffusion algorithm (Bruenn 1985) in the multiple ID 
"ray-by-ray" transport approximation As do Marek 
& Janka (2009), the ORNL team employs the lower-K, 
softer nuclear EOS. They have found robust explosions 
for a variety of progenitors (Bruenn et al. 2007,2010; 
Mezzacappa et al. 2007) and have suggested that the 
mechanism is a combination of neutrino heating, nuclear 
burning (of infalling oxygen at the shock), the SASI, and 
inelastic neutrino-electron and neutrino-nucleon scatter- 
ing. However, the latter should be a subdominant 5—10% 
effect. The potential role of nuclear burning on infall had 
been studied and eliminated earlier (Burrows & Lattimer 
1985; H.-T. Janka, private communication) and though 
included in many other previous simulations has not been 
seen to have played a positive role. Infalling oxygen 
should burn long before reaching the shock wave that 
has stalled near ~150— 200 kilometers and, moreover, 
the amount of nuclear energy available is rather small. 
Why the results of these two groups differ so substan- 
tially (when their approaches are ostensibly so similar), 

^ Also used by the Marek &; Janka (2009), and first introduced 
into supernova theory in more primitive form by Burrows, Hayes, 
& Fryxell (1995). 



in both explosion energy when explosions are claimed, 
and in whether they see explosions at all, is a puzzle. 

Those using the fully 2D radiation/hydrodynamics 
code VULCAN (Livne 1993; Livne et al. 2004,2007; Bur- 
rows et al 2006, 2007ab), apart from obtaining neutrino- 
driven explosions in the contexts of a 8.8-M0-progenitor 
(Burrows, Dessart, & Livne 2007) and accretion-induced- 
collapse (Dessart et al. 2006b), do not witness neutrino- 
driven explosions for most progenitors. VULCAN/2D 
incorporates an implicit, time-dependent, multi-energy- 
group transport scheme that has both multi- angle (S^) 
and flux-limited (MGFLD) variants in two spatial dimen- 
sions. VULCAN is the flrst and (to date) only multi- 
group supernova code to operate in 2D in both the hy- 
dro and radiation sectors. Interestingly, using a multi- 
angle variant of VULCAN, Ott et al. (2008) have shown 
that for non-rotating models the results for the multi- 
angle and MGFLD simulations are similar. VULCAN 
also incorporates "2.5-D" MHD and, hence, is the only 
code currently employed in supernova studies with both 
multi-group transport and MHD capabilities. With it. 
Burrows et al. (2007c) have explored MHD-driven jet 
explosions (LeBlanc & Wilson 1970; Bisnovatyi-Kogan, 
Popov, & Samokhin 1976; Akiyama & Wheeler 2005; 
Wheeler & Akiyama 2007), but concluded that such ex- 
plosions, which require very rapid rotation at bounce 
that is very unlikely to be generic, might obtain only 
in the rare "hypernova" case and could not explain the 
typical supernova. 

However, using VULCAN and after simulating for 
around a second after bounce. Burrows et al. (2006, 
2007a) observed vigorous dipolar g-mode oscillations 
that damp by the asymmetric emission of sound waves 
that steepen into shock waves. This acoustic power has 
been adequate to explode all progenitors, but would be 
aborted if the neutrino mechanism, or some other mech- 
anism, obtained earlier. This acoustic mechanism is con- 
troverial, takes a long time (many seconds) to achieve ex- 
plosion energies of ~10^^ ergs, and has not been seen by 
other groups. However, no other group has calculated for 
the physical time necessary to witness vigorous g-mode 
excitation. Unlike all VULCAN simulations, most other 
simulations have been done with the inner core in lE|fl 
As a result, we suggest that other groups have not sim- 
ulated long enough to see this phenomenon. Moreover, 
we suspect that installing a spherical inner core will par- 
tially suppress the excitation of a core g-mode that is fun- 
damentally dipolar, not spherical, and whose amplitude 
should be largest at the center. A solid criticism of the 
acoustic mechanism comes from Weinberg & Quataert 
(2008), who using approximate analytic techniques, sug- 
gest that the non-linear excitation of dissipative daughter 
modes via a parametric resonance might render the am- 
plitudes of these dipolar g-mode oscillations too small to 
power a supernova. Unfortunately, the wavelengths of 
these daughter modes are smaller than reasonable grid 
spacings and can not easily be captured by extant su- 
pernova codes. Hence, whether a very delayed acoustic 
power mechanism for explosion is at all viable in any 
circumstance remains to be seen. 

However, the delayed neutrino mechanism seems com- 

^ The VULCAN modelers used an enabling Cartesian grid in the 
inner 30—50 kilometers. 
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pelling, if only it can be shown to work robustly 
and to give canonical explosion energies of ~10^^ ergs 
generically. Burrows & Goshy (1993) suggested that 
the neutrino-driven explosion is a critical phenomenon. 
There is a critical luminosity at a given mass accretion 
rate into the shock above which there is no accretion 
solution and the mantle explodes. Recently, Murphy & 
Burrows (2008) have shown that this critical luminos- 
ity decreases and the neutrino-matter coupling efficiency 
increases by ~30% in going from ID to 2D. Murphy & 
Burrows (2008) have diagnosed the increase in heating 
efficiency in 2D (vis a vis ID) as the increase in the aver- 
age dwell time (Thompson, Quataert, & Burrows 2005) 
of matter in the gain region due to neutrino-driven con- 
vection and the SASI, before it settles into the cooling 
region and onto the inner core. What hasn't been shown, 
and is our goal with this paper, is the corresponding ef- 
fect in going from 2D to 3D. 

The absence in VULCAN of inelastic scattering terms 
and any corrections for general relativistic effects might 
compromise its conclusions and be the reason those us- 
ing VULCAN are not seeing neutrino-driven explosions 
in the generic case, or when other groups obtain such 
explosions (however tepid). However, the inelastic terms 
are small, and relativity's effects nearly cancel. A more 
plausible explanation for the ambiguity rife in the field 
and the marginality (and rarity) of neutrino-driven ex- 
plosions among published simulations is that the neu- 
trino mechanism may not work well in 2D. It may be that 
the neutrino mechanism is truly viable only in 3D and 
that 3D effects have been the missing ingredient needed 
to explain supernova explosions with delayed neutrino 
driving. This is the thesis of our paper. We posit that 
the viability of the neutrino-driven mechanism for core- 
collapse supernova explosions and the energy of explosion 
(when they explode) are monotonic functions of dimen- 
sion. We note that it is only recently, with the advent 
of computers of sufficent size and speed, that theorists 
have been able to test such an hypothesis. Furthermore, 
we suggest that the ambiguity of the extant 2D results 
is a symptom of its marginality in 2D. Moreover, as did 
Murphy & Burrows (2008), we conjecture that the in- 
creased dwell time in the gain region in going to 3D, due 
both to the extra degree of freedom for non-radial mo- 
tion and the character of the cascade in 3D turbulence 
to small turbule size^ renders the heating efficiency and 
all the other ancillary quantities that support explosion 
adequate to transform "marginality and ambiguity" into 
"robustness and viability" (qil). 

3. CASTRO CODE AND GENERAL METHODOLOGY 

To address these issues in a full 3D context, we 
have employed CASTRO (Compressible ASTRO; Alm- 
gren et al. 2010), a new multi-dimensional radia- 
tion/hydrodynamic code. CASTRO is an Eulerian, 
structured grid, compressible, radiation/hydrodynamics 
code that incorporates adaptive mesh refinement (AMR). 
It uses a second-order unsplit piecewise-parabolic (PPM) 
Godunov methodology for general convex equations of 
state. The most general treatment of self-gravity in CAS- 
TRO uses multigrid to solve the Poisson equation for the 

^ In 2D, the turbulent cascade is artifically inverted and flows to 
large scales. 



gravitational potential. For the calculations presented 
here, the monopole approximation is used. 

AMR in CASTRO uses a nested hierarchy of rectan- 
gular grids that are simultaneously refined in both space 
and time. CASTRO uses a recursive integration proce- 
dure in which coarse grids are advanced in time, fine grids 
are advanced multiple steps to reach the same time as 
the coarse grids, and the data at different levels are then 
synchronized. For regridding, an error estimation proce- 
dure evaluates where additional refinement is needed and 
grid generation procedures dynamically create or destroy 
fine grid patches to achieve the desired local resolution. 
CASTRO is implemented within the BoxLib framework 
that handles data distribution, communication, memory 
management, and I/O for parallel architectures. 

To most closely parallel the approach of Murphy & 
Burrows (2008), while generalizing to three dimensions, 
our calculations were done with CASTRO using the sim- 
ple neutrino heating and cooling algorithm described in 
Murphy & Burrows (2008). We solve the fully compress- 
ible equations of hydrodynamics: 

dtp=-V-{pu) (1) 
dt (pu) = -V • (puu) - + pg 
dt {pE) = -V • {puE) +pu + pu • g + p(H - C) , 

where p, T, p, g, and u are the fluid density, temper- 
ature, pressure, gravitational acceleration, and velocity. 
The total energy is given hy E = e -\- ^u'^ ^ where e repre- 
sents the internal energy. The equation of state provides 
closure and is a function of p, T, and electron fraction, 
Yg. While a sophisticated nuclear equation of state is 
used (Shen et al. 1998), the neutrino heating and cooling 
rates behind the stalled shock wave via the super-allowed 
charged- current reactions involving free nucleons assume 
a given electron neutrino and anti-electron neutrino lumi- 
nosity (taken to be the same). This luminosity is varied 
from simulation to simulation, but is held constant dur- 
ing a simulation. The neutrino heating, and cooling, 
C, rates are those derived in Janka 2001, used in Murphy 
& Burrows (2008), and given by 

n = lM4xlO^^ (---^ tVt^^Vx (2) 

yiQ52 erg s-1 J \4 MeV J ^ ^ 

\ r J ^ [gs_ 

and 

(3) 

where Lj^^ is the electron neutrino driving luminosity, 
Tiy^ is the electron neutrino temperature, r is the distance 
from the center of the star, and Yp are the neutron and 
proton fractions and Tj^^ is the electron neutrino optical 
depth. Note that these approximations assume that the 

This approach, using eqns. ([2| & (|3| in place of full 
transport, enabled the extensive (ana computationally 
expensive) parameter study in ID, 2D, and 3D we report 
here. The single progenitor we employ for all our runs is 
the non-rotating, solar-metallicity 15-Mq red-supergiant 
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model of Woosley & Weaver (1995), but we need only 
the inner ~ 10,000 km for our simulations of core col- 
lapse. The calculations commence at the onset of infall 
and are carried out to and after bounce. The explo- 
sion time is formally determined when the average shock 
radius reaches 400 km and does not recede during sub- 
sequent evolution (see Table 1). Many of the ID and 2D 
simulations are carried to beyond a second after bounce 
(up to ~1.4 seconds). The 3D simulations were all car- 
ried out to beyond 400 milliseconds after bounce. The 
Liebendorfer et al. (2005) scheme for determining the 
electron fraction on infall was used, which involves tying 
Ye to the mass density achieved for densities below the 
trapping density. 

Interior to a radius of 200 km, all of our simulations 
have zones smaller than a kilometer. Exterior to 200 km, 
the resolution dynamically follows the shock via adaptive 
grids. In 3D, our Cartesian domain is a cube of length 
10,000 km. The domain is discretized at the coarsest level 
with a 304^ uniform grid, and several levels of adaptive 
refinement are used so that cells at the finest level are 
a factor of 64 finer in each direction. This results in a 
resolution of ~0.5 km (which we refer to as the "effec- 
tive" resolution) in the interior 200 km of the star and 
in other regions at the highest refinement level. For our 
axisymmetric simulations, we employ a 2D domain that 
is 10,000 km by 5,000 km, discretized with a uniform 
coarse grid of 256 by 128 cells. All of our 2D simulations 
cover the full 180° angular range from north to south 
pole. Employing the same refinement criteria used in 3D 
yields an effective resolution of ~0.6 km for our 2D simu- 
lations. For our spherically symmetric (ID) simulations, 
we can obtain, and have tested, higher resolution runs. 
Additional resolution leads to minimal differences (^10 
ms) in the reported explosion times. The ID results pre- 
sented in this paper are for a radial domain of 256 coarse 
cells (with similar refinement to that employed in the 
multi-D simulations), distributed over 5,000 km for an 
effective resolution of ~0.3 km. 

4. RESULTS 

Depicted in Fig. [l]are critical curves in driving lumi- 
nosity versus accretion rate (M) through the shock for 
calculations performed in ID, 2D, and 3D, all else be- 
ing equal. The luminosity is in units of 10^^ ergs s~^ 
and the mass accretion rate is in solar masses per second 
(Mq s~^). M is 4:7Tr'^ p\vii\ in the infalling material just 
exterior to the shock wave. Above each curve the core 
explodes and below each curve it does not. As shown in 
Murphy & Burrows (2008), the ID curves approximately 
recapitulate the ID analytic theory found in Burrows & 
Goshy (1993). The important result in this paper is the 
position of the 3D curve. It is ~15— 25% below the corre- 
sponding 2D curve, which is in turn ~30% below the ID 
curve. Moreover, the advantage of going to 3D is larger 
for higher driving lumniosities. The upshot is that the 
magnitude of the driving neutrino luminosity necessary 
to "supernova" a given core for a given mass accretion 
rate through the shock is reduced in going from ID to 3D 
by ~40— 50%, a rather large, perhaps enabling, effect. 

Note that in our study each run is performed for a given 
neutrino and anti-neutrino luminosity, while the accre- 
tion rate onto the inner core evolves with time, running 



from high to low values of M in a way determined by the 
initial inner density profile of the chosen progenitor star 
(in this case the 15-M0 progenitor of Woosley & Weaver 

1995). Hence, the M is not fixed during a run and if 
another progenitor model had been chosen the tempo- 
ral evolution of M would have been different. When the 
M reaches the value for a given fixed Lj^^ at which the 
core explodes, we identify this Lj^^—M pair as a point 
on the critical curve (Fig. [T]). Performing this exercise 
for a range of luminosities maps out the corresponding 
critical curve in whatever number of dimensions is being 
studied. 

Since before explosion the trajectory in Liy^{t)—M{t) 
space followed by a given progenitor is a weak function 
of dimension, depending mostly on the initial progen- 
itor density profile, which in turn determines M and 
thereby the inner core density, temperature, and Yg pro- 
files (given an EOS, a neutrino transport algorithm, and 
quasi- hydrostatic equilibrium), changing the dimension 
of a simulation is not met with significant compensating 
shifts in the actual Ly^(t) vs. M(t) trajectory with which 
to compare our critical curves. The result is an undimin- 
ished and monotonic advantage at higher dimensions of 
significant magnitude. The effect of dimension is seen to 
be not merely a few percent, but nearly a factor of two 
in a central aspect of the collapse context, the driving 
neutrino luminosity. This new result leads us to suggest 
that lack of access to 3D computational capabilities has 
been a major retarding factor in progress during the last 
few decades towards the solution to the core-collapse su- 
pernova problem. 

Plotted in Fig. [2 are various curves depicting the tem- 
poral evolution of the average radius of the shock (in kilo- 
meters) for various driving electron neutrino luminosities 
and for simulations in ID, 2D, and 3D. This figure shows 
that for luminosities for which the ID and 2D model 
shocks remain stagnant for long periods, the 3D models 
explode much earlier. For instance, at Ly^ = 1.9 x 10^^ 
ergs s~^, in ID the core doesn't explode even after 1.4 
seconds, the 2D core explodes after around ~0.8— 1.0 sec- 
onds, but the 3D core explodes after only ~250 millisec- 
onds. In all cases, the time to explosion (if there is an 
explosion) is shorter at higher dimension than at lower 
dimension. Bruenn et al. (2007,2010) witness explo- 
sions in all their recent 2D simulations, but these results 
have not been reproduced by others (cf. Marek & Janka 
2009). Moreover, all their explosions pile up at similar 
early times, as does the 3D simulation they are currently 
tending. This suggests that whatever is causing their ex- 
plosions does not much distinguish between 2D and 3D 
in the way we so clearly see in our suite of simulations. 
The reason for this is unclear, but we note that when 
we obtain early explosions (in this paper, at the high- 
est driving neutrino luminosities), the difference in the 
time to explosion in 2D and 3D is similarly significantly 
reduced. 

Table 1 clearly demonstrates that the time to explosion 
is shorter at higher dimension than at lower dimension 
and provides a more extended compilation of ID, 2D, and 
3D exploding models and the approximate times at which 
they explode. The table is arranged so that overlapping 
horizontal rows, though done for a different number of 
dimensions, have the same driving luminosities. This 
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format clearly reveals that, all else being equal, the time 
to explosion is significantly shorter at higher dimension. 
For example, the ID model at L^^ = 2.5 x 10^^ ergs s~^ 
explodes around ^0.75 seconds while the corresponding 
2D model explodes near 0.2 seconds. Interestingly, this 
is the time at which the 3D model at the much lower 
luminosity of 1.9 x 10^^ ergs s~^ explodes. The time it 
takes the 1.9 x 10^^ ergs s~^ model in 3D to lift the av- 
erage shock radius from '^200-300 kilometers to ~1200 
kilometers during the early explosion phase is ~200 mil- 
liseconds, at which point the shock is moving at a speed 
of ~30,000 km s~^. As indicated on the figure, early in 
the incipient explosion phase the average shock radius 
gradually, but steadily, accelerates. 

In the panels in Fig. |3]we compare representative en- 
tropy color maps of simulations for the same neutrino lu- 
minosities and times after bounce, but for different num- 
bers of dimensions. The top two panels contrast ID (left) 
and 2D (right) runs, both for a luminosity of 2.5 x 10^^ 
ergs s~^ and at ~468 milliseconds after bounce. Note 
that while the 2D run is exploding, the ID run is not, 
though the physical model and inputs are otherwise iden- 
tical. The positions of the shocks in these simulations, as 
indicated by the abrupt color transition, are very differ- 
ent. The bottom panels compare models in 2D (left) and 
3D (right), both for a luminosity of 1.9 x 10^^ ergs s~^ 
and at ~422 milliseconds after bounce. While the shock 
in the 2D simulation is still stalled, the shock in the 3D 
simulation has launched dramatically. These snapshots 
together and in conjunction demonstrate pictorially the 
qualitative dependence on dimension of the outcome of 
collapse. Since supernovae occur in three spatial dimen- 
sions, we conclude that this dimensional dependence is of 
direct relevance to the viability of the neutrino heating 
mechanism for core-collapse supernova explosions. More- 
over, the magnitude of the "advantage" of going to 3D 
(see Fig. [T| should dwarf that of any refinements in the 
microphysics, such as inelastic neutrino-matter scatter- 
ing, nuclear burning upon infall, or general relativity. 
The latter are at most "~10%" effects and are small com- 
pared with the ~40-50% effects we identify here that are 
associated with dimension, in particular in going to 3D. 

Figure |4] depicts entropy scatter plots near and af- 
ter the time of explosion of the corresponding higher- 
dimensional run for 1) ID and 2D models (top) and 2) 
2D and 3D models (bottom). In each case, the lower- 
dimensional model of a given comparison (either top or 
bottom) has not exploded during the time spanned by 
the plots. The top panels are for L^^^ = 2.1 x 10^^ ergs 
s"-*^ and the bottom panels are for hy^ = 1.9 x 10^^ ergs 
s~^. The times to explosion are given in Table 1. Figure 
|4] clearly indicates that, all else being equal, higher en- 
tropies are achievable at higher dimension. Figure [5] de- 
picts the corresponding mass-weighted average entropy 
in the gain region versus time for the same ID (black), 
2D (blue), and 3D (red) models and makes the same 
point. Importantly, the average entropy for the 3D run 
in the runup to explosion is ~1.5 units higher than in 
2D. This expands the physical extent of the gain region, 
makes it less bound, and creates what seems to be a more 
"explosive" situation (cf. Figs, ij [2j and[3|. The higher 
average entropies can be traceatonigher average dwell 
times of a Lagrangian parcel of matter in the gain region 



(Murphy & Burrows 2008). We plan to present a more 
detailed analysis of this and other effects in a subsequent 
paper. 

In addition, as Figs. [3] and [6] indicate, even though 
we do see pre-explosion and blast asymmetries in 3D, 
we do not see a strong i = 1 axially-symmetric "SASI" 
oscillation, such as is obtained in 2D in this paper and 
by others (Blondin et al. 2003; Bruenn et al. 2007,2010; 
Burrows et al. 2006,2007a; Buras et al. 2006ab; Marek & 
Janka 2009). Rather, the free energy available to power 
instabilites seems to be shared by more and more de- 
grees of freedom as the dimension increases (Iwakami et 
al. 2008). In ID, as demonstrated by Murphy & Burrows 
(2008), one sees a strong radial oscillation when the lumi- 
nosity is near critical. However, 2D models do not mani- 
fest this radial mode, but instead execute an £ = 1, m = 
dipolar oscillation. In 3D, this dipolar oscillation exists 
(Blondin & Mezzacappa 2007; Fernandez 2010), but in 
the linear limit competes with the m = { — 1, 1} modes. 
The ^ = l,m = mode is even less in evidence in the 
non-linear limit, which is achieved early in our simula- 
tions. Hence, the strong dipolar symmetry seen in 2D 
and identified as the fundamental characteristic of the 
"SASI" may not survive in 3D. At the very least, the 
^ = 1, m = mode does not dominate in our non-rotating 
3D models. Whether rotation can change this conclusion 
remains to be seen (Guilet, Sato, & Foglizzo 2010). 

5. CONCLUSIONS 

In this paper, we have performed ID, 2D, and 3D nu- 
merical simulations of the collapse, bounce, and explo- 
sion of a massive star core to isolate the effect of spatial 
dimension in the context of the neutrino heating mech- 
anism of core-collapse supernova explosions. We have 
found that both the viability of explosion (measured by 
the driving neutrino luminosity needed to overcome a 
given mass accretion rate at the shock and its associ- 
ated accretion pressure tamp) and the position of the 
"critical curve" are monotonically increasing functions 
of dimension, with 3D more viable than 2D by 15—25% 
and 3D more viable than ID by almost a factor of two. 
Some had thought that 3D would prove similar to ID, but 
our calculations do not support this expectation. More- 
over, we have discovered that the time to explosion is 
significantly shorter in 3D than 2D, all else being equal. 
These results suggest that a key missing ingredient in the 
recipe for this astrophysically central phenomenon has 
been access to 3D simulation tools and that the neutrino- 
driven explosion mechanism is fundamentally 3D. Two- 
dimensional simulations have to date yielded marginal 
explosions and/or ambiguous answers (for a discussion 
see Janka et al. 2007). 

Furthermore, we find that the prominent £ = 1, m = 
dipolar mode of shock oscillation, often identified in 2D 
as a central aspect of the "SASI," is little in evidence 
in 3D. The free energy available to excite hydrodynamic 
instabilities seems instead to be shared by more degrees 
of freedom, diminishing the amplitude of this axisym- 
metric dipolar component. Whether this conclusion is a 
function of progenitor and/or rotation has yet to be de- 
termined, but the preliminary indications are that this 
sloshing mode, so visible in 2D simulations, does not sur- 
vive as a central feature of core collapse when three spa- 
tial dimensions are allowed for. 
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We emphasize that in this study, in order to isolate 
the crucial effect of dimension, we employed a very sim- 
ple neutrino "transport" approach. This has enabled 
us to discover what we suggest is a key aspect of the 
core-collapse supernova phenomenon. However, detailed 
neutrino transport must be incorporated to determine 
the true systematics with stellar progenitor of the ex- 
plosion energies, nucleosynthesis, residual neutron star 
(and black hole) masses, pulsar kicks, and blast mor- 
phologies. If we had included such transport at this 
stage in the development of computer hardware and the 
computational arts, a single 3D simulation would have 
required many, many years of continuous execution on 
the largest existing and available massively-parallel plat- 
forms. As a result, the insights we have achieved via 
our more modest study might have been delayed or ob- 
scured for years. This highlights the continuing benefits 
of a multi-pronged, varied, and flexible strategy to ad- 
dress this most important, though refractory, astrophys- 
ical problem. 
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TABLE 1 

Explosion times and accretion rates. The explosion time, texp, is calculated when the average shock radius, (i^shock), 

REACHES 400 KM AND DOES NOT RECEDE DURING SUBSEQUENT EVOLUTION. THE ACCRETION RATE AT THE TIME OF EXPLOSION IS 
CALCULATED JUST EXTERIOR TO THE SHOCK WHERE THE INFALLING ENVELOPE IS SPHERICAL, AND IS GIVEN BY Mexp = 4:7Tpvr'^. A " — " 
INDICATES THE MODEL DID NOT EXPLODE DURING THE ALLOTTED MAXIMUM SIMULATION TIME (1.4 SECONDS). THE BOLDED TIMES ARE 
THOSE FOR MODELS WHICH BOTH EXPLODED WITHIN 1.4 SECONDS AND WERE CALCULATED IN EITHER BOTH ID AND 2D OR BOTH 2D AND 
3D. As THIS TABLE INDICATES, THERE WERE MODELS THAT WERE CALCULATED IN ID, 2D, AND 3D, BUT THE CORRESPONDING ID MODELS 

DID NOT EXPLODE WITHIN 1.4 SECONDS. 
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Fig. 1. — Critical curves in electron-neutrino driving luminosity (Lj^^) versus accretion rate (M) for calculations performed in ID (black), 
2D (blue), and 3D (red). The luminosity is in units of 10^"^ ergs s"-*^ and the mass accretion rate (M) is in solar masses per second. 
Importantly, we include results in 3D. Note that the driving electron neutrino luminosity is always accompanied by an associated anti- 
electron neutrino luminosity (Lj^g)- See ^for a discussion of the meaning of this plot and its salient features. 
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Fig. 2. — Depicted here are curves showing the temporal evolution of the average radius of the shock (in kilometers) for various driving 
electron neutrino luminosities and for simulations in ID, 2D, and 3D. The ID models are rendered as black dashed lines, the 2D models as 
solid blue lines, and the 3D models as thick solid red lines. Each line is indexed by the corresponding electron-neutrino luminosity of the 
simulation, in units of 10^^ ergs s"-*^. Implied is an equal anti-electron-neutrino luminosity. Time is in seconds after bounce. See the text 
in ^for a discussion of this figure. 
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Time = 0.468 s 




Fig. 3. — In these panels we compare representative entropy maps of two simulations for the same driving neutrino luminosities and times 
after bounce, but for different numbers of dimensions. The top two panels contrast ID (left) and 2D (right) runs, both for an electron-type 
neutrino luminosity of 2.1 x 10^^ ergs s"-*^ and at ~0.468 milliseconds after bounce. The bottom panels compare models in 2D (left) and 3D 
(right), both for an electron- type neutrino luminosity of 1.9 x 10^-^ ergs s"-*^ and at 0.422 milliseconds after bounce. The same colormap is 
used for all four panels. Note that while the two top panels are for the same luminosity and epoch after bounce, only the 2D simulation has 
exploded. Similarly, while the two bottom panels are at the same luminosity and time after bounce, the development of the 3D simulation 
is qualitatively different from that of the corresponding 2D run. Note also the different general morphologies of the 3D (bottom) and 2D 
(top) models that explode. See the text in ^and ^for a discussion. 
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Fig. 4. — Scatter-plot entropy profiles versus radius for two different driving luminosities, comparing ID and 2D (top) and 2D and 3D 
(bottom). The ID (red)-2D (black) comparison at the top is for L^^g = 2.1 x 10^^ ergs and the 2D (red)-3D (black) comparison at the 
bottom is for Ljy^ = 1.9 x 10^-^ ergs s"-*^. For both top and bottom two different times after bounce are given — the times for both left 
panels are near the onset of the explosion of the associated higher-dimensional model. Note that the peak entropies and average entropies 
are always higher for the higher-dimensional runs, and that for a given dimension the higher the driving luminosity the higher the peak 
entropy achieved. Note also that different horizontal radius scales are employed for the left and the right panels and that the associated 
model names are given on the plots. See the discussion in ^ 
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Fig. 5. — The mass-weighted average entropy in the gain region versus time after bounce (in seconds) for the ID(black), 2D(blue), 
and 3D(red) models with L^^^ = 1.9 x lO^"^ ergs s"-*^. This figure, representative of corresponding figures at other driving luminosities, 
demonstrates the higher entropies achieved in the gain region behind the shock in going to 3D. After ~0.2 seconds after bounce the 
difference between the 2D and 3D model-average entropies is ~1.5 units. See ^for a discussion of the meaning and relevance of this figure. 
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Fig. 6. — Blast Morphology. These panels depict an evolutionary sequence (top left, top right, bottom left, bottom right) of the 
Li/g = 1.9 X 10^^ erg s"-*^ model in 3D. The top left panel depicts a phase well before the onset of explosion. The two surfaces in each panel 
are p = 10^^ g cm""^ (blue interior) and Yg = 0.47 (outer) for four different times after bounce (0.156, 0.201, 0.289, and 0.422 seconds). The 
scale is more than 2000 kilometers on a side. Note that the crude axis of the explosion is not along any of the three Cartestian directions 
and that there is no obvious i = \ asymmetry. See ^for a discussion of the possible implications of this figure and of Fig. [3] for the 
relevance of the = 1, m = 1 "SASI" mode in the context of core-collapse supernovae. 



